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Abstract. 

Using density-functional molecular dynamics simulations we analyzed the cooling-rate 
effects on the physical properties of GeS2 chalcogenide glasses. Liquid samples were cooled 
linearly in time according to T(t) = Tg —yt where y is the cooling rate. We found that our 
model leads to a promising description of the glass transition temperature T g as a function of y 
and gives a correct T g for experimental cooling rates. We also investigated the dependence of 
the structural properties on the cooling rate. We show that, globally, the properties determined 
from our simulations are in good agreement with experimental values and this even for the 
highest cooling rates. In particular, our results confirm that, in the range of cooling rates 
studied here, homopolar bonds and extended charged regions are always present in the glassy 
phase. Nevertheless in order to reproduce the experimental intermediate range order of the 
glass, a maximum cooling rate should not be exceeded in numerical simulations. 



PACS numbers: 61.43.Fs,61.43.Bn,71.23.-k,71.15.Pd 
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1. Introduction 

During the last fifteen years the unceasingly growing interest for chalcogenide glasses has 
lead to numerous works. This is in particular true for germanium disulfide glasses for which a 
large amount of experimental U |2] [3] |4) and theoretical J5] |6] [7] studies have been carried 
out. Indeed its known properties used for example in optical amplifiers, memory switching 
devices or anti-reflection coatings [9] make him a good candidate for intensive research. 
Among the different research means, molecular dynamics (MD) simulations can be a very 
interesting tool to provide detailed information on the physical properties of such glassy 
systems. Firstly because they allow to investigate the structure in full microscopic detail 
giving access to the position of the atoms and, secondly because they are useful to study 
dynamical phenomena accessible to such simulations e.g., for time scales between 10~ 13 and 
1(T 8 s. 

If the temperature of a liquid is decreased so much that the relaxation time of the system 
becomes longer than the time scale of the computer simulation or of the experiment, the 
system undergoes a kinetic arrest and, provided that it does not crystallize, will undergo a glass 
transition and remain trapped in a disordered configuration. However it has been demonstrated 
in both, experiment IflOl [TTI and computer simulations Ifl2l [TJl , that the properties of the 
resulting glass, like the density or the glass transition temperature, will depend on its thermal 
history and in particular on the rate at which the sample is cooled down. 
Previous studies ||5] |6] UJ have already validated our "Cook and Quench" model to produce 
and study GeS2 glasses using approximate ab initio molecular dynamics simulations. 
Nevertheless, because of the time scale of our computer simulations which is many orders 
of magnitude shorter than the typical experimental one, the glass transition temperature, for 
example, appears to be significantly higher than the one observed in the laboratory (^750K 
fl2j EJUS)) and thus it is necessary to see how the properties of the so-obtained glass depend 
on the way it was produced. 

Thus, in the present paper, we focus on the cooling rate effects on some physical properties of 
glassy GeS2- Firstly we investigate how it does affect the glass transition temperature T g and 
secondly we study how it does affect the structural properties of the glassy samples. 
The paper is consequently organized as follows. In section|2]we briefly present the theoretical 
model used in our calculations, results and discussions are presented in section [3] and finally 
in section |4] we summarize the major conclusions of our work. 

2. Theoretical framework 

Computations were performed using Fireball96, an approximate ab initio molecular dynamics 
code based on the local-orbital electronic structure method developed by Sankey and 
Niklewski lfl6l . The electronic structure is described using density functional theory (DFT) 
ifTTl within the local density approximation (LDA) lfl8l and the non-local pseudo-potential 
scheme of Bachelet, Hamann and Schliiter |fl9l . To reduce the CPU time we used the non- 
selfconsistent Harris functional ll20l with a set of four atomic orbitals ( 1 "s" and 3 "p") per 
atom that vanish outside a cutoff radius of 5ao (2.645A). This model has been successfully 
used the last ten years for several different chalcogenide systems ll5l l6l 171 12T1 l22l l23l . 
All the calculations of the present simulations were performed in the microcanonical 
ensemble, with a time step At=2.5fs and using only the F point to sample the Brillouin 
zone. The initial configuration of our system was a crystal of oc-GeS2 in a cubic cell of 
19.2lA containing 258 atoms with standard periodic boundary conditions, melted at 2000 
K for 60ps in order to obtain an equilibrium liquid. Then, we relaxed this system further 
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at 2000K for 50ps and we choose five different liquid samples (approximately every lOps) 
during this process. Each sample was then quenched down to 300K through the glass 
transition temperature T g . Quenches were carried out by a linear velocity rescaling according 
to T(t) = Tq — yr where y is the cooling rate. Six different cooling rates were used: 3, 0.6, 0.3, 
0.12, 0.09 and 0.06x 10 14 K/s. At 300K each sample was relaxed over 50ps, i.e., 20000 time 
steps. Configurations were saved every 20 steps and the results were averaged for each sample 
over these 1000 configurations. Furthermore we averaged the results for the 5 samples of each 
cooling rate thus all the data presented below have been averaged over 5 x 1000 configurations, 
and the error bars (when given) represent the usual standard deviation. 

3. Results 

3.1. Glass transition 

As usually done in simulations we tried to identify the glass transition of our samples by 
representing the evolution of the average potential energy as a function of the temperature 
during the quench (an example is proposed for one of the samples cooled at 0.09 x 10 14 K/s 
(Fig. [TJ), the glass transition being localized by the change in the evolution of the energy. The 
glass transition temperature T g is determined by the intersection of the two linear regressions 
at high and low temperature (Fig. [T). For each cooling rate we averaged the values of the T g 
obtained for each of the five GeS2 samples, thus we obtained an evolution of the average T g 
with the cooling rate. As shown in Fig. |2]the error bars remain huge indicating an insufficient 
sampling (it is worth noting that the computer time needed to perform the simulations for the 
lowest values of y is of the order of 12 weeks for one sample) and it is therefore difficult to 
extract an accurate description of the evolution of T g with the cooling rate. Nevertheless we 
attempted to fit these average simulation values assuming a power-law dependence of T g with 
the cooling rate as suggested by the mode-coupling theory ll24ll : 

T g (y)=T c +(Ay)^ (1) 

with T c = 676 ±75K, A = 3.0 ±2.3 x 10 31 and 8 = 17.1 ±0.1 (the errors given for the 
parameters of the fit have been evaluated by fitting the extreme (lowest and highest) Tg values 
enclosed in the error bars). We observed that a variation of y by about 1 decade gives rise 
to a variation of T g of about 75K which is not that larger, compared to the difference of 
magnitude of the cooling rates, than the variation of 10K measured in real experiments for 
different materials flOl [TTl l24l . Then by extrapolating the results of our fit down to usual 
experimental cooling rates i.e. 10°-10 5 K/s, we observed at this scale that a variation of yby 
about 1 decade gives rise to a variation of T g of about 13K, which is in agreement with the 
experimental variation of 10K previously mentioned. We have represented in Fig. |2]the result 
of this work together with a few experimental data l2l [T4l \l5l . The agreement between the 
extrapolation of the fit and the experimental values of T g shows that our model is able to give 
a correct tendency of the variation of T g with y, even though the poor statistics prevents us 
from having accurate estimates of T g at a specific (high) cooling rate. This can be improved 
with more simulations in order to reduce the error bars. 

3.2. Structural properties 

3.2.1. Radial pair correlation functions and bonding properties 

In glassy GeS2 the basic building blocks are GeS4 tetrahedra, connected together forming 
a random network. The structural disorder is reflected by the absence of long range order 
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and by the wide distribution of bond lengths and bond angles. Structural information may be 
extracted from the radial pair correlation function g(r) which can be defined for a given a, (3 
pair by: 

^> = i5^ § 8 <'-'<>> (2) 

where p is the number density of the system, c a the fraction of species a in the system, i the 
atoms of species a and j the atoms of species p. For each cooling rate we averaged the radial 
pair correlation functions g a p(r) of the five samples, so we are able to compare the evolution 
of the average gap(r) according to the cooling rate (Fig. [3landl4li. 

The bond lengths appear to be in good agreement with experimental data(2| since we find 
2.23A for the Ge-S bond (expe.: 2.2lA), 2.9lA and 3.49A for respectively the edge and 
corner sharing Ge-Ge connections (expe: 2.91 A and 3.42A). This good agreement is true even 
for the highest cooling rates. The main influence of the cooling rate is reflected in the small 
peak corresponding to homopolar bonds between 2.2-2.6A for Ge (Fig. [3j> and 2.1-2.45A 
for S (Fig. @). Our results indicate that the number of homopolar bonds decreases with the 
cooling rate and it is therefore justified to address the question of the existence of homopolar 
bonds at experimental cooling rates which is still an open question (Cai and Boolchand using 
Raman scattering experiment, found the existence of homopolar bonds in glassy GeS2 l25ll 
while Petri and Salmon found no evidence of such bonds in gGeS2 using neutron diffraction 
lfT4l studies). According to Fig. [3] and Fig. |4]it seems that the decrease of the proportion 
of homopolar bonds slows down for the lowest cooling rates and tends towards a limit of 
respectively 1.9% for the Ge atoms and 1.2% for the S atoms. These limiting values are small 
but nonzero and therefore our simulation results seem to confirm the existence of homopolar 
bonds in experimental glassy GeS2- 

The simulation gives access to the positions of the atoms, therefore we can also obtain 
information on the connectivity of the network. In our approach we focused on the 
ratio between edge and corner sharing tetrahedra and the evolution of the different local 
environments of the Ge (Tab. [TJ and S atoms (Tab. 



y(10 14 K/s) 


Proportion of Ge atoms [%] 
Ge(S 4 ) Ge(S 3 ) Ge(GeS 3 ) 


3 


93.0 ±2.0 


2.3 ±1.8 


3.7 ±2.1 


0.6 


95.0 ±1.6 


1.7 ±1.1 


1.9 ±2.0 


0.3 


96.2 ±1.8 


0.3 ±0.3 


2.8 ±2.0 


0.12 


96.4 ±1.8 


1.7 ±1.8 


1.9 ±2.0 


0.09 


97.6 ±1.7 


0.1 ±0.1 


1.8 ±1.0 


0.06 


98.1 ±1.0 


0.0 ±0.0 


1.9 ±1.1 



Table 1. Evolution of the local structural environment of Ge atoms as a function of the cooling 
rate y 

First it is worth noting that the proportion of edge and corner-sharing links is almost a 
constant independent of the cooling rate: 84%±1.8 of corner-sharing and 16%±1.8 of edge- 
sharing bounds, values in good agreement with experimental data (2]. This is in contrast with 
the evolution of the proportions of Germanium (TabQ} and Sulfur (Tabf2]i in their standard 
environment (respectively a four fold S coordination for Ge and a two fold Ge coordination 
for S) which appreciably increase with decreasing cooling rate. The second point concerns, 
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y(10 14 K/s) 


S(Ge 2 ) 


Proportion of S atoms[%] 
S(Ge) S(Ge 3 ) S(GeS) 


S(Ge 2 S) 


3 


67.9 ±2.2 


14.3 ±1.5 


12.9 ±1.8 


3.1 ±1.3 


1.8 ±0.9 


0.6 


71.0 ±4.1 


12.9 ±2.0 


12.6 ±1.8 


2.1 ±0.8 


0.4 ±0.3 


0.3 


75.5 ±2.6 


11.4 ±1.0 


10.7 ±1.2 


2.0 ±1.3 


0.4 ±0.6 


0.12 


76.1 ±2.7 


10.9 ±1.2 


10.5 ±1.3 


1.7 ±0.8 


0.6 ±0.6 


0.09 


78.2 ±3.3 


10.4 ±1.6 


9.9 ±1.4 


1.1 ±1.1 


0.4 ±0.3 


0.06 


78.5 ±1.0 


10.4 ±0.4 


9.8 ±1.0 


0.8 ±0.3 


0.4 ±0.3 



Table 2. Evolution of the local structural environment of S atoms as a function of the cooling 
rate y 



as expected, the decrease of the chemical disorder with decreasing cooling rate. Indeed the 
proportion of under-coordinated Ge atoms (2.3% for the fastest cooling rate) disappears for 
the slowest cooling rate. And the proportions of non-bridging S atoms and over-coordinated 
S atoms, respectively 14.25% and 12.85% decrease to 10.4% and 9.8%. 
These results indicate that the cooling rate has an impact on the structure of the glass. 
Nevertheless while certain types of structural "defects" disappear at low rate, others survive 
and can therefore be considered as inherent of the glassy structure. 

3.2.2. Neutron static structure factor 

An alternative way to analyze the structure is to compute the static neutron structure fac- 
tor S(q) which can be directly compared to neutron scattering experiments: 

m = ^b j b k (e i ^-r^ (3 ) 

where N is the number of atoms, and bj is the neutron scattering factor for atom j. 
As for the radial pair distribution functions, we averaged for each cooling rate the total neu- 
tron structure factors of the five samples, permitting thus the comparison of the evolution of 
the average structure factors as a function of the cooling rate (Fig. |5J. 

First we note the accurate description of glassy GeS2 reflected in the good agreement between 
the simulated curves and the experimental one, and this even for the highest cooling rates. Al- 
though the simulated and experimental curves present differences in the range 0-2.5 A" 1 it has 
already been shown J6l|2l|26) that the physical properties of the so-simulated glassy samples 
are in good agreement with experiment. The First Sharp Diffraction Peak (FSDP), signature 
of the intermediate range order (IRO) in amorphous samples, appears at «lA~' and is, as 
expected in such simulations, slightly underestimated J5J. Size effects can be considered as 
an explanation, nevertheless we note that a decrease of the cooling rate globally improves the 
calculated structure factor and in particular the FSDP (Fig. [5]), which highlights that the IRO 
is also cooling rate dependent. However there is no linear/regular evolution of S(q) with the 
cooling rate. To illustrate the improvement of S(q) in the FSDP region with the decrease of the 
cooling rate, we represent (Fig. [6]) the difference, between and 2.5A -1 , of the experimental 
total neutron structure factor and: 

1) the total neutron structure factor of the liquid state at 2000K: 

S(q)ex P -S(q)n q = AS(q) exp -i iq 
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2) the average total neutron structure factor of the simulated glass quenched down to 300K 
at the rate y: 

S(q)exp -S(q) y = AS(q) exp - y 

The AS(q) exp -iiq is a reference representing the biggest variation between the experimental 
and simulated S(q) and should be compared to AS(q) exp -y. We see in Fig. |6]that by decreasing 
the cooling rate, the difference AS(q) exp -y decreases. In particular it appears that for cooling 
rates higher than 0.3xl0 14 K/s the differences AS(q) exp -y are very close to the difference 
AS(q)exp-Hq ■ F° r the highest cooling rate (3xl0 14 K/s) this difference is even quasi identical 
to AS(q) exp -Uq which indicates that there is no real change in the IRO between the fastest 
quenched glass and the liquid phase. One can thus argue that the fastest cooled samples are 
too similar to the liquid and can therefore not be considered as glassy GeS2 samples. This 
defines a limit for the maximum cooling rate usable in MD simulations in order to avoid 
interferences between the liquid and the glassy state. In our simulations this limit appears to 
be between 0.6xl0 14 K/s and 0.3xl0 14 K/s. Nevertheless this limit is directly related to the 
method i.e. ab-initio MD simulations and indirectly related to our model i.e. the nature of the 
glass and the characteristics of the atomic pseudopotentials and therefore the numerical value 
of this limit can not be straightforwardly extended to other simulated glassy systems. 

3.2.3. Atomic charges 

Even if atomic charges cannot actually be determined experimentally, relevant tools such as 
Lowdin 11271 or Mulliken 11281 population analysis can be used to compare different configura- 
tions with the same description. In the present work the Lowdin description has been chosen 
in order to compare the dependence of the atomic charges on the cooling rate. It should how- 
ever be mentioned that the non-self-consistent Harris functional is known to overestimate the 
charge transfers between the atoms. 

The atomic charge q is calculated by the difference between the number of electrons of the 
neutral atom and the "real" number of electrons of the atom in the glass. We found no depen- 
dence of the Lowdin charges on the cooling rate. Thus it is necessary to correlate the charges 
with the evolution of the proportion of each atomic type in its local structural environment 
with the cooling rate (Tab. Q]and|2]i. 

As expected the general polarity of the Ge-S bond is found with a charge transfer in an or- 
dered Ge(S4)i/2 configuration of +0.94 for the Ge atoms and -0.46 for the S atoms. As it has 
already been shown in a previous work [6,| Ge charges are always positive and decrease with 
the number of neighbors whereas S charges are more variable with respect to the local envi- 
ronment: from strongly negative charges for non-bridging S atoms (-1.07) to almost neutral 
charges for 3 fold Ge-coordinated S atoms. And even if the existence of positively charged S 
atoms in the environment q| >S-S is confirmed in our present work, Tab. [2] shows that this 
kind of local structural environment disappears rapidly with decreasing cooling rate. These 
structures obtained at high cooling rate may be explained by the results shown in Tab. [3] 
representing the charge of an atom in a given local structural environment and the proportion 
of atoms in this environment in an equilibrium GeS2 liquid at 2000K. Indeed it appears that 
the configuration of the liquid is quite similar, at least for the S atoms, to the fastest cooled 
glass configurations (Tab. This result indicates that the fastest cooled glasses are, in the 
literal sense, frozen liquids. This confirms what we have already detected in the total neutron 
structure factor. 

Positively and negatively charged zones inside the glass have been reported in our previous 
study [6 1 for the highest cooling rate. In order to see if the atomic charges measured for the 
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Environment q Proportion 

Ge(S 4 ) +0.96±0.04 64.83% 

Ge(S 3 ) +1.02±0.06 27.41% 

Ge(GeS 3 ) +0.75±0.1 1.14% 

Ge-S-Ge -0.43±0.15 60.22% 

Ge-S -0.92±0.21 26.03% 

ll >S-Ge -0.03±0.11 9.74% 

Ge-S-S -0.1 ±0.1 8 2.07% 

>S-S 0.32±0.2 1.60% 



Table 3. Lowdin charges according to the local structural environments in an equilibrium 
liquid GeS 2 at 2000K 



lowest cooling rate confirm or reject the existence of such zones we looked at the short-range 
charge deviation AQsr of a particle i defined by: 

n(<) a( i) 

A& S (i)=?(i) + ll (4) 
P <J) 

This allows to take into account the atomic charge q(i) of a given particle ; as well as the 
charges on its n(i) nearest neighbors (determined from the radial pair distribution function). 
Whereas for a crystalline structure, in which no bond defects are present, this value is almost 
zero for all the particles, positive and negative values appear for glassy samples. 
In Tab. [4] we have reported the evolution of the number of charged zones (a charged zone 
contains at least 2 nearest neighbors having the same sign for AQsr) as a function of the 
cooling rate. 

The number of positively charged zones (particles with AQsr > +0.3) is almost constant and 



y(10 14 K/s) 


Averaged number of charged zones 
AQsr > +0.3 AQ SR < -0.3 


3 


9.0 ±1.0 


20.8 ±1.2 


0.6 


9.0 ±2.0 


19.2 ±0.8 


0.3 


8.4 ±1.6 


17.2 ±0.8 


0.12 


8.6 ±1.4 


16.0 ±1.0 


0.09 


7.8 ±1.8 


15.8 ±0.2 


0.06 


9.2 ±0.8 


16.4 ±0.6 



Table 4. Evolution of the average number of charged zones as a function of the cooling rate y 



equal to 9, with no dependence on the cooling rate, whereas the number of negatively charged 
zones (particles with AQsr < -0.3) decreases slightly from ~20 for the samples cooled at 
3 x 10 14 K/s to ~16 for the samples cooled at a rate smaller than 0.3 x 10 14 K/s (this is again 
coherent with the idea of a maximum cooling rate usable in MD simulations). The negatively 
charged zones are principally made of Ge atoms coordinated to one or more non-bridging 
sulfur atoms, and for the fastest cooled glasses, of a few S-S-Ge structures. As already shown 
(Tab. [2]) these structures disappear with decreasing cooling rate. This observation correlated 
to the diminution of the proportion of non-bridging sulfur atoms gives an explanation to the 
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14 Averaged number of atoms per charged zones 
Y(10 K/s) AQsR>+03 AQsR < _Q.3 

3 5.7 ±0.1 2.1 ±0.0 

0.6 5.2 ±0.4 2.1 ±0.1 

0.3 5.2 ±1.1 2.0 ±0.0 
0.12 5.2 ±0.2 2.1 ±0.0 
0.09 5.2 ±1.7 2.1 ±0.0 
0.06 4.2 ±0.2 2.1 ±0.0 

Table 5. Evolution of the average number of atoms per charged zones as a function of the 
cooling rate y 



decrease of the number of negatively charged zones with the cooling rate for rates higher than 
0.3 x 10 14 K/s. The positively charged zones are exclusively made of Ge atoms linked to over- 
coordinated S atoms. 

In addition we have reported the average number of atoms per charged zone as a function 
of the cooling rate (Tab. [5]). This shows that the size of the positively and the negatively 
charged zones (respectively with 5 and 2 atoms per zone) is independent of the cooling rate 
and therefore remains constant. It is worth noting that the global neutrality of the glass is 
always respected. Our results show that the existence of charged zones in glassy GeSi is not 
influenced by the variation of the cooling rate. They confirm thus those of our previous study, 
and show that extended charged zones (whose manifestation has also been detected recently 
for other chalcogenide systems by Taraskin et al. 11291 ) reflect the broken chemical order of 
the glass and are therefore inherent to the amorphous state. 

4. Conclusion 

Through DFT based molecular dynamics simulations we have analyzed the effect of the 
cooling rate on some properties of glassy GeS2- Influence of the cooling rate on the glass 
transition temperature as well as on the structural properties has been studied. 
Due to a lack of statistics mainly due to computer time limitations (especially for the lowest 
cooling rates), the detailed variation of T g with the cooling rate could not be obtained. 
Nevertheless the extrapolation of our results to "realistic" cooling rates is in good agreement 
with the experimental glass transition temperature. 

Analyzing the radial pair distribution functions and the local structural environments for each 
cooling rate, we find that the number of S and Ge homopolar bonds as well as the number 
of coordination defects decrease with the cooling rate. Nevertheless the decrease of the 
homopolar bonds seems limited and therefore it is reasonable to think that this type of defect is 
present in real glasses. However calculations at lower cooling rates should be done to confirm 
this observation. 

The study of the simulated total neutron structure factor has confirmed the reliability of our 
model in the description of glassy GeS2- We have analyzed the effect of the cooling rate on 
the intermediate range order whose signature is the first sharp diffraction peak. The simulated 
FSDP is closer to the experimental one for the slowest cooled glasses. By comparing with 
the properties of the liquid state, we have shown that a maximum cooling rate should not be 
exceeded in the simulation in order to reproduce the IRO characteristic of the glassy phase. 
The value of this maximum cooling rate will depend on the details of the model used to 
describe a given system. The existence of a maximum cooling rate has also been supported 
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by the analysis of the charges which has revealed that the electronic configuration of the 
fastest cooled glasses is close to the one obtained in liquid GeS2- In addition the existence of 
positively and negatively charged regions in the amorphous state has been clearly confirmed 
even for the lowest cooled samples and seems therefore inherent to the glassy state. These 
regions will have an important impact on the properties of samples containing metallic ions 
as shown recently l26l . 
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Figure 2. Glass transition temperature as a function of the cooling rate, fit and extrapolation 
to very low cooling rates - comparison with experiment [1 I Ti l 1 51 . 




Figure 3. Partial radial pair distribution function g(r)[Ge-Ge] as a function of the cooling 
rate. Enlargement corresponds to details of the peak due to homopolar bonds. 




Figure 4. Partial radial pair distribution function g(r)[S-S] as a function of the cooling rate. 
Enlargement corresponds to details of the peak due to homopolar bonds. 




Figure 5. Total simulated static neutron structure factors. Enlargement corresponds to the 
FSDP part of the S(q). 
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Figure 6. AS(q) for the IRO, between experiment and liquid GeS2, and between experiment 
and glassy GeS2 for the different cooling rates. 



